! Characteristics Calculation
subroutine CalcChars()
use Data_Mesh, only: nCells, meshFaceOwner, &
                     nBoundaryFields, boundaryMeshName, boundaryMeshStart, boundaryMeshSize, &
                     meshFaceCenters
use DataCells, only: CharsCellPresent
use DataFaces, only: CharsFacePresent, fluxesFacePresent, FaceVarsNew, FaceVarsOld, Acoef, CharsNew
use Data_Cfoam, only: btype, utmp, Ufon, RoFon,Pfon,Tfon, ViscPartOnly
use DataPerfectGas, only: DensityPT
use DataTime, only: TimeCell, dtCell
implicit none

    integer :: iCell, iBoundField, iFace, iCellOwner
    logical :: CalcAcoef
    real(8) :: NoSlipCoef
    real(8) :: U0(3),P0,Temp,Pres,Rho,Energy
    real(8) :: R(3), Q,dQdT,c, x, TimeX, Time, rMod
    real(8) :: y,z, Umean, Lr, Px

    do iCell = 1,nCells

        if(CharsCellPresent(iCell)) cycle

        ! The calculation of chars in cell iCell
        call CalcCellChars(iCell)

        CharsCellPresent(iCell) = .true.
    end do

    ! Calculation charachteristics for boundary faces
    do iBoundField = 1,nBoundaryFields
        if (boundaryMeshName(iBoundField)(1:4).eq.'proc') cycle

        if (btype(iBoundField).eq.1) then ! Outlet
            do iFace=1+boundaryMeshStart(iBoundField),boundaryMeshStart(iBoundField)+boundaryMeshSize(iBoundField)
                if(fluxesFacePresent(1,iFace)) cycle

                Acoef(2,iFace) = Acoef(1,iFace)
                CharsNew(1:5,2,iFace) = CharsNew(1:5,1,iFace)

                cycle

!                if(CharsCellPresent(iCellOwner)) cycle

                FaceVarsNew(5,2,iFace) = utmp(5,iBoundField)
                Temp = utmp(4,iBoundField)
                FaceVarsNew(2,2,iFace) = UFON
                FaceVarsNew(3:4,2,iFace) = 0d0
                FaceVarsNew(1,2,iFace) = DensityPT(FaceVarsNew(5,2,iFace),Temp)
                
                FaceVarsNew(1:5,2,iFace) = FaceVarsNew(1:5,1,iFace)

                CalcAcoef = .false.
                call CalcInOutletBoundaryChars(iFace,CalcAcoef)
            end do
        end if

        if (btype(iBoundField).eq.2) then ! Inlet
            do iFace=1+boundaryMeshStart(iBoundField),boundaryMeshStart(iBoundField)+boundaryMeshSize(iBoundField)
                if(fluxesFacePresent(1,iFace)) cycle

                if(ViscPartOnly)then
                    x = meshFaceCenters(1,iFace)
                    y = meshFaceCenters(2,iFace) !- 0.5d0
                    z = meshFaceCenters(3,iFace)

                    Umean = 0.01
                    Lr = 0.005
                    Px = - 12d0 * 1.5e-5 / Lr**2 * Umean

                    utmp(5,iBoundField) = Pfon + Px * x
                    UFON = -0.5d0 / 1.5e-5 * Px * z*(Lr-z)
                    utmp(4,iBoundField) = Tfon
                end if

                FaceVarsNew(5,2,iFace) = utmp(5,iBoundField)
                Temp = utmp(4,iBoundField)
                FaceVarsNew(2,2,iFace) = UFON
                FaceVarsNew(3:4,2,iFace) = 0d0
                FaceVarsNew(1,2,iFace) = DensityPT(FaceVarsNew(5,2,iFace),Temp)

                FaceVarsOld(1:5,2,iFace) = FaceVarsNew(1:5,2,iFace)

! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Boundary conditions from subroutine CalcQ3D
!                iCellOwner = meshFaceOwner(iFace)
!                r(1:3) = meshFaceCenters(1:3,iFace)
!                rMod = sqrt(sum(r(1:3)*r(1:3)))
!                c = sqrt(1.4*PFon/RoFon)
!                time = TimeCell(iCellOwner)+dtCell(iCellOwner)
!                call CalcQ3D(r,time,Q,dQdt,p0,U0)
!!                p0 = PFon
!!                U0(1) = 10d0
!!                U0(2:3) = 0d0
!
!                Pres = p0
!                Rho = RoFon + (P0-PFon)/c**2
!                Energy = P0 / (0.4 * (RoFon + (P0-PFon)/c**2)) + 0.5d0 * sum(U0(1:3)*U0(1:3))
!
!                FaceVarsNew(1,2,iFace) = Rho
!                FaceVarsNew(2:4,2,iFace) = U0(1:3)
!                FaceVarsNew(5,2,iFace) = Pres
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


                if(sum(U0(1:3)*r(1:3)).gt.0)then
                    CalcAcoef = .true.
                else
                    CalcAcoef = .false.
                end if
                call CalcInOutletBoundaryChars(iFace,CalcAcoef)
            end do
        end if

        if (btype(iBoundField).eq.3) then ! Symetry plane ! wall, no-Slip
            do iFace = 1+boundaryMeshStart(iBoundField),boundaryMeshStart(iBoundField)+boundaryMeshSize(iBoundField)
                if(fluxesFacePresent(1,iFace)) cycle

                NoSlipCoef = 0d0
                call CalcWallBoundaryChars(iFace,NoSlipCoef)
            end do
        end if

        if (btype(iBoundField).eq.4) then ! Symetry plane ! wall, full-Slip
            do iFace = 1+boundaryMeshStart(iBoundField),boundaryMeshStart(iBoundField)+boundaryMeshSize(iBoundField)
                if(fluxesFacePresent(1,iFace)) cycle

                NoSlipCoef = 1d0
                call CalcWallBoundaryChars(iFace,NoSlipCoef)
            end do
        end if
    end do

end subroutine CalcChars

! Calculation of charachteristics for inlet or outlet
! important:
! FaceVarsNew(1:5,2,iFace) - boundary values
subroutine CalcInOutletBoundaryChars(iFace,CalcAcoef)
use Data_Mesh, only: meshSf
use DataFaces, only: FaceVarsNew, Acoef, CharsNew
use DataPerfectGas, only: ACoefPRho
implicit none
    integer,intent(in) :: iFace
    logical,intent(in) :: CalcAcoef

    real(8) :: n(3),normal(3)
    real(8) :: a(3,3)
    real(8) :: Pres,Rho,U(3),Ut(3)
    real(8) :: Chars3,Sound

    n(1:3) = meshSf(1:3,iFace)
    normal(1:3) = n(1:3)/dsqrt(sum(n(1:3)*n(1:3)))
    ! transformation {x,y,z} -> {normal,tan1,tan2}
    call coor(normal(1),normal(2),normal(3),a)
!    a(1:3,1:3) = 0d0
!    a(1,1) = 1d0
!    a(2,2) = 1d0
!    a(3,3) = 1d0
!    call transform(FaceVarsNew(2:4,2,iFace),Ut(1:3),a)
     
    Pres = FaceVarsNew(5,2,iFace)
    Rho = FaceVarsNew(1,2,iFace)
    U(1:3) = FaceVarsNew(2:4,2,iFace)

    if(CalcAcoef) then
        ACoef(2,iFace) = ACoefPRho(Pres,Rho) ! GammY * (Pres * Rho**(-Gamma))**GammPower
    else
        Acoef(2,iFace) = Acoef(1,iFace)
    end if

    call Vars2Chars(U,Pres,Rho,iFace,a,Acoef(2,iFace),Chars3,Sound,CharsNew(1:5,2,iFace))

end subroutine CalcInOutletBoundaryChars

! calculation of charachteristics for wall
subroutine CalcWallBoundaryChars(iFace,NoSlipCoef)
use Data_Mesh, only: meshFaceOwner
use DataFaces, only: CharsNew, Acoef, FaceVarsNew
use DataPerfectGas, only: EntropyPRho
implicit none
    integer,intent(in) :: iFace
    real(8),intent(in) :: NoSlipCoef

    integer :: iCellOwner
    real(8) :: Pres,Rho

    iCellOwner = meshFaceOwner(iFace)
    CharsNew(1,1,iFace) = -CharsNew(2,1,iFace)
    CharsNew(4:5,1,iFace) = NoSlipCoef * CharsNew(4:5,1,iFace)

    Acoef(2,iFace) = Acoef(1,iFace)

    CharsNew(1:2,2,iFace) = CharsNew(1:2,1,iFace)
    CharsNew(4:5,2,iFace) = CharsNew(4:5,1,iFace)

    Pres = FaceVarsNew(5,1,iFace)
    Rho = FaceVarsNew(1,1,iFace)
    CharsNew(3,2,iFace) = EntropyPRho(Pres,Rho)
    
    FaceVarsNew(1:5,2,iFace) = FaceVarsNew(1:5,1,iFace)

end subroutine CalcWallBoundaryChars

!Definitioin of charachteristics for cell iCell
subroutine CalcCellChars(iCell)
use Data_Mesh, only: meshCells
use DataFaces, only: fluxesFacePresent
implicit none
    integer,intent(in) :: iCell

    integer :: iFaceNum, iFace

    do iFaceNum = 1,6
        iFace = meshCells(iFaceNum,iCell)

!        if(fluxesFacePresent(1,iFace).or.fluxesFacePresent(2,iFace)) cycle

        call CalcFaceCellChars(iCell,iFace)
    end do

end subroutine CalcCellChars

!Definitioin of charachteristics in cell iCell for face iFace
subroutine CalcFaceCellChars(iCell,iFace)
use Data_Mesh, only: meshFaceOwner, opposingOwnerFaceLabel, opposingNeighbourFaceLabel, &
                     meshSf, dF, dB
use DataCells, only: SubStepCellVars
use DataFaces, only: ACoef, CharsOld, CharsFacePresent
use DataPerfectGas, only: PressureAllVars, AcoefPRho
use Data_Cfoam, only: epsilonSmoothed, epsilonPanikovski
implicit none
    integer,intent(in) :: iCell,iFace

    integer :: iCellOwner,iCellNeighbour, oppFace, iSide, oppSide, iChar
    real(8) :: AcoefF, AcoefB
    
    real(8) :: EigenVal3Opp,SoundOpp,CharsOpp(5)
    real(8) :: EigenVal3SC,SoundSC,CharsSC(5)
    real(8) :: EigenVal3C,SoundC,CharsC(5), Epsilon
    real(8) :: EigenVal3Old,SoundOld
    real(8) :: EigenVal1SC, EigenVal2SC
    
    real(8) :: ds
    
    real(8) :: n(3),normal(3),a(3,3)
    real(8) :: U(3),Pres,Rho,Energy,Temp

    iCellOwner = meshFaceOwner(iFace)
    if(iCellOwner.eq.iCell)then
        oppFace = opposingOwnerFaceLabel(iFace)
        iSide = 1
        ds = -dB(iFace)
    else
        oppFace = opposingNeighbourFaceLabel(iFace)
        iSide = 2
        ds = dF(iFace)
    end if

!    if(CharsFacePresent(iSide,iFace)) return

    ! Definition of Characteristics, Eigen Values and Sound Velocities in
    Pres = PressureAllVars(SubStepCellVars(1:5,iCell))
    Rho = SubStepCellVars(1,iCell)
    ACoef(iSide,iFace) = ACoefPRho(Pres,Rho) ! GammY * (Pres * Rho**(-Gamma))**GammPower

    n(1:3) = meshSf(1:3,iFace)
    normal(1:3) = n(1:3)/dsqrt(sum(n(1:3)*n(1:3)))
    ! transformation {x,y,z} -> {normal,tan1,tan2}
    call coor(normal(1),normal(2),normal(3),a)

    ! Opposite Face
    if(iCell.eq.meshFaceOwner(oppFace))then
        oppSide = 1
    else
        oppSide = 2
    end if
    call FaceVars2Chars(oppFace,iFace,a,Acoef(iSide,iFace), EigenVal3Opp,SoundOpp,oppSide,CharsOpp(1:5))
    ! Cell Center, on the Time SubStep
    call SubStepCellVars2Chars(iCell,iFace,a,Acoef(iSide,iFace), EigenVal3SC,SoundSC,CharsSC(1:5))
    ! Cell Center, on the Old Time Step
    call CellVars2Chars(iCell,iFace,a,Acoef(iSide,iFace), EigenVal3C,SoundC,CharsC(1:5))

    call FaceVars2Chars(iFace,iFace,a,Acoef(iSide,iFace), EigenVal3Old,SoundOld,iSide,CharsOld(1:5,iSide,iFace))

    ! Definition of Panikovskii-Viscosity:
    epsilon = epsilonSmoothed
    if(iSide.eq.1)then
        ! RIGHT CELL
        if(EigenVal3Old.le.0 .and. EigenVal3Opp.gt.0)epsilon = epsilonPanikovski !0.7 ! empirical calculation methods for important points ! SK
        if(EigenVal3Old+SoundOld.le.0 .and. EigenVal3Opp+SoundOpp.gt.0)epsilon = epsilonPanikovski !=0.7
        if(EigenVal3Old-SoundOld.le.0 .and. EigenVal3Opp-SoundOpp.gt.0)epsilon = epsilonPanikovski !=0.7
    else
        ! LEFT CELL
        if(EigenVal3Old.ge.0 .and. EigenVal3Opp.lt.0)epsilon = epsilonPanikovski !0.7 ! empirical calculation methods for important points ! SK
        if(EigenVal3Old+SoundOld.gt.0 .and. EigenVal3Opp+SoundOpp.le.0)epsilon = epsilonPanikovski !=0.7
        if(EigenVal3Old-SoundOld.gt.0 .and. EigenVal3Opp-SoundOpp.le.0)epsilon = epsilonPanikovski !=0.7
    end if

    EigenVal1SC = EigenVal3SC + SoundSC
    EigenVal2SC = EigenVal3SC - SoundSC

    ! Definition of Characteristics on new time level
    ! CharsNew(1:5,iSide,iFace)
    call CalcChar(EigenVal1SC,2)
    call CalcChar(EigenVal2SC,1)

    do iChar = 3,5
        call CalcChar(EigenVal3SC,iChar)
    end do
    
!    CharsFacePresent(iSide,iFace) = .true.

contains

    ! Defenition of Characteristics on New Time Level, Extrapolation of Characteristics
    subroutine CalcChar(EigenVal,i)
    use DataTime, only: dtCell,dt2Cell
    use DataFaces, only: CharsNew
    implicit none
        real(8),intent(in) :: EigenVal
        integer,intent(in) :: i

        real(8) :: Qrhs, Qmax,Qmin, Qnew

        ! Definition of Inhomogeneous Right Hand Side
        Qrhs = (CharsSC(i) - CharsC(i)) / dt2Cell(iCell) + &
               EigenVal * (CharsOpp(i) - CharsOld(i,iSide,iFace)) / ds

        ! Min and Max Definition
	    Qmax = dmax1(CharsOld(i,iSide,iFace),CharsSC(i),CharsOpp(i)) + dtCell(iCell) * Qrhs
	    Qmin = dmin1(CharsOld(i,iSide,iFace),CharsSC(i),CharsOpp(i)) + dtCell(iCell) * Qrhs

        ! CABARET equation, extrapolation
        Qnew = (2d0 * CharsSC(i) - CharsOpp(i)*(1.0d0-epsilon)) / (1.0d0+epsilon)

        ! Conservative Correction, Maximum Principle
        Qnew = min(Qnew,Qmax)
        Qnew = max(Qnew,Qmin)

        CharsNew(i,iSide,iFace) = Qnew

    end subroutine CalcChar

end subroutine CalcFaceCellChars

! Analytic 1D solution
subroutine CalcQ(r,time,Q,dQdt,p,U)
use Data_cfoam, only: PFon,RoFon,UFON
implicit none
    real(8),intent(in) :: r(3),time
    real(8),intent(out) :: Q,dQdt, p,U(3)

    real(8) :: Pi,rMod,PiR4,c,ek,a, RC

    Pi = 3.1415926535897932384626433832795
    rMod = 0d0 ! sqrt(sum(r(1:3)*r(1:3)))
    PiR4 = 1d0 ! 4.0d0 * Pi * rMod
    c = sqrt(1.4*PFon/RoFon)
    RC = r(1) / c
    ek = 1d0 !*0.1d0
    a = 40d0 !c / 100.d0 ! 
!    a = 40.0d0
    Q = ek * sin(2.0d0 * a * Pi * time) / PiR4
    dQdt = ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * time) / PiR4

!    Q = ek * sin(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
!    dQdt = ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4

    Q = 100d0 * max(0d0,min(time,0.01d0)) * ek * sin(2.0d0 * a * Pi * max(0d0,time - RC)) / PiR4
    dQdt = 100d0 * max(0d0,min(time,0.01d0)) * ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * max(0d0,time - RC)) / PiR4

    Q = -ek * sin(2.0d0 * a * Pi * (time - RC)) / PiR4
    dQdt = -ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * (time - RC)) / PiR4

!    Q = 0d0
!    dQdT = 0d0
    
    P = dQdt + Pfon
    U(1) = dQdt / (c*RoFon) + UFON ! + Q / rMod ) / (rMod * ROFon)

!    P = Pfon
!    U(1) = 50d0 ! -max(0d0,min(time-0.01d0,0.01d0))
!    U(2:3) = 0d0

end subroutine CalcQ

! Analytic 3D solution
subroutine CalcQ3D(r,time,Q,dQdt,p,U)
use Data_cfoam, only: PFon,RoFon,Ufon
implicit none
    real(8),intent(in) :: r(3),time
    real(8),intent(out) :: Q,dQdt, p,U(3)

    real(8) :: Pi,rMod,PiR4,c,ek,a

    Pi = 3.1415926535897932384626433832795
    rMod = sqrt(sum(r(1:3)*r(1:3)))
    PiR4 = 4.0d0 * Pi * rMod
    c = sqrt(1.4*PFon/RoFon)
    ek = 0.1d0
    a = 5.0d0
!    Q = ek * sin(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
!    dQdt = ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4

    Q = ek * sin(2.0d0 * a * Pi * (time - rMod / c)) / PiR4
    dQdt = ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * (time - rMod / c)) / PiR4
!
!    Q = 100d0 * max(0d0,min(time-0.01d0,0.01d0)) * ek * sin(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
!    dQdt = 100d0 * max(0d0,min(time-0.01d0,0.01d0)) * ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
!
!    Q = 100d0 * max(0d0,min(time,0.1d0)) * ek * cos(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
!    dQdt = -100d0 * max(0d0,min(time,0.1d0)) * ek * 2.0d0 * a * Pi * sin(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
!
!    Q = ek * cos(2.0d0 * a * Pi * time) / PiR4
!    dQdt = - ek * 2.0d0 * a * Pi * sin(2.0d0 * a * Pi * time) / PiR4

!    if(time.lt.1d-10)then
!        Q = 0d0
!        dQdt = 0d0
!    end if

    P = dQdt + Pfon
    U(1:3) = r(1:3) * (dQdt / c + Q / rMod) / (rMod * ROFon) !  + 1d0 / (rMod * PiR4)
    
!    return
!    
!    P = Pfon
!    U(1) = Ufon
!    U(2:3) = 0d0
!    
!!    U(1:3) = r(1:3) * (Q / rMod) / (rMod * ROFon) !  + 1d0 / (rMod * PiR4)
!
!!    P = Pfon
!!    U(1:3) = 10d0 * max(0d0,min(time,0.1d0)) / 0.1d0 * r(1:3) / rMod !* (dQdt / c + Q / rMod) / (rMod * ROFon) !  + 1d0 / (rMod * PiR4)
!
!!    U(1) = 7d0
!!    U(2:3) = 0d0
!
!    call CalcQ(r,time,Q,dQdt,p,U)

end subroutine CalcQ3D

! Analytic 2D solution
subroutine CalcQ2D(r,time,Q,dQdt,p,U)
use Data_cfoam, only: PFon,RoFon
implicit none
    real(8),intent(in) :: r(3),time
    real(8),intent(out) :: Q,dQdt, p,U(3)

    real(8) :: Pi,rMod,PiR4,c,ek,a

    Pi = 3.1415926535897932384626433832795
    rMod = sqrt(sum(r(1:3)*r(1:3)))
    PiR4 = 4.0d0 * Pi * rMod
    c = sqrt(1.4*PFon/RoFon)
    ek = 10d0
    a = 10.0d0
    Q = ek * sin(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
    dQdt = ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4

    Q = ek * sin(2.0d0 * a * Pi * (time - rMod / c)) / PiR4
    dQdt = ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * (time - rMod / c)) / PiR4

!    Q = 100d0 * max(0d0,min(time-0.01d0,0.01d0)) * ek * sin(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
!    dQdt = 100d0 * max(0d0,min(time-0.01d0,0.01d0)) * ek * 2.0d0 * a * Pi * cos(2.0d0 * a * Pi * max(0d0,time - rMod / c)) / PiR4
    
    P = dQdt + Pfon
    U(1:3) = r(1:3) * (1d0 + dQdt / c + Q / rMod) / (rMod * ROFon) !  + 1d0 / (rMod * PiR4)

!    P = Pfon
!    U(1:3) = 500d0 * max(0d0,min(time-0.01d0,0.01d0)) * r(1:3) !* (dQdt / c + Q / rMod) / (rMod * ROFon) !  + 1d0 / (rMod * PiR4)

end subroutine CalcQ2D

! Analytic 2D solution
subroutine CalcQ2D_Hankel(r_in,time_in,p,U,rho)
!use Data_cfoam, only: PFon,RoFon, UFON
!USE IFPORT, only: DBesY0, DBesJ0, DBesJn,DBesYn
!implicit none
!    real(8),intent(in) :: r_in(3),time_in
!    real(8),intent(out) :: p,U(3), rho
!
!    integer :: i
!    real(8) :: A, omega, c, k, s, omegaT, dr, rMod, rMod1, drMod !, r1(3)
!    real(8) :: BesY0s, BesJ0s, r(3)
!!    real(8) :: gradBesY0(3),gradBesJ0(3)
!    real(8) :: gradBesY0(3),gradBesJ0(3), BesJ1s, BesY1s, ddBesY0(3),ddBesJ0(3)
!    real(8) :: Ut(3), gradP(3), Err(3), Pt, divU
!    real(8) :: gradBesY01, gradBesJ01 !, ddBesY0, ddBesJ0
!    real(8) :: s1, ErrPt,ErrPtt,Ptt,laplaceP, ErrPt1
!    real(8) :: Pxx,Pyy,Pzz, dY1,dJ1,dY1c,dJ1c, s2, dJ0s,dY0s, dJ0s1,dY0s1
!    real(8) :: time, M, beta
!    real(8) :: Phi, PhiS,PhiT,PhiX,PhiY
!    
!    real(8) :: beta2, Mx1c
!    
!    real(8) :: Tau,dTauDt,dTauDx, dTauDy, D
!    real(8) :: r1,r2, Rn, dRnDt, dRndX, dRndY
!    real(8) :: J0,Y0, dJ0Dr,dY0Dr, cosT,sinT
!    real(8) :: Mr, dMrDt, dMrDx, dMrDy, Mr1
!    
!    real(8) :: HQ, dHdrQ, HdQdTau
!    real(8) :: dJ0ds,dY0ds
!    
!    real(8) :: dr1dX, dr1dy, dr1dT    
!    
!    real(8) :: r1_,r2_,Rn_,Mr_,dMrDtau_
!
!    r = r_in
!    time = time_in
!    r(3) = 0d0
!
!    A = 5d0
!    omega = 1d0
!    c = sqrt(1.4*PFon/RoFon) ! c = 340d0
!    !r(1) = r(1) + Ufon * time
!    r(1) = r(1) + Ufon * time
!    k = omega/c
!    M = UFON / c
!    rMod = sqrt(sum(r*r))
!    beta2 = 1-M**2
!
!    Mx1c = M*r(1)/c
!
!    D = (r(1) - Ufon*time)**2 + beta2*r(2)**2
!    if(D.lt.0)then
!        write(*,*) '@@Fatal Error: D<0'
!        Stop '@@Fatal Error: D<0'
!    end if
!    D = sqrt(D)
!    
!    Tau = (c*time - M*r(1) - D) / (c*beta2)
!     
!!!    dTauDt = (1d0 - (time*(1-beta2) + Mx1c) / D) / beta2
!!!    dTauDx = (M/c - (M/c*time + r(1)/c**2) / D) / beta2
!!!    dTauDy = - r(2) / (c**2 * D)
!!
!!    dTauDt = (1d0 - (time*M**2 - Mx1c) / D) / beta2
!!!    dTauDx = (M/c - (-M/c*time + (beta2+1d0) * r(1)/c**2) / D) / beta2
!!    dTauDx = (-M/c - (r(1)/c**2-M/c*time) / D) / beta2
!!!    dTauDx = -Mr / (Ufon*(1d0-Mr))
!!    dTauDy = - r(2) / (c**2 * D)
!
!
!!    r1 = r(1)+Ufon*Tau
!    r1 = r(1)-Ufon*Tau
!    r2 = r(2)
!    Rn = c*(time-Tau) ! sqrt(r1**2 + r2**2) !
!
!    Mr = M * r1 / Rn
!
!!    Mr1 = 1-Mr
!!    Mr1 = 1+Mr
!    Mr1 = 1-Mr
!    
!
!    dTauDt = 1d0 / Mr1 ! (1d0 - (time*M**2 - Mx1c) / D) / beta2
!    dTauDx = Mr / (Ufon * Mr1) ! (-M/c - (r(1)/c**2-M/c*time) / D) / beta2
!    dTauDy = r2 / (c*Rn * Mr1) ! - r(2) / (c**2 * D)    
!!    dRndX = r1 / Rn
!!    dRndY = r2 / Rn
!
!!    dRndX = -c * dTauDx
!!    dRndY = -c * dTauDy
!!    dRndTau = -c ! Ufon * dRndX
!    
!    dRndX = -c * dTauDx ! r1 * (1d0 - Ufon*dTauDx) / Rn ! 
!    dRndY = -c * dTauDy ! (r2 - r1*Ufon*dTauDy) / Rn !
!    dRnDt = c*(1-dTauDt) ! -r1 * Ufon / Rn ! -c ! Ufon * dRndX
!    
!    dr1dX = 1d0 - Ufon * DtauDx
!    dr1dy = -Ufon * DtauDy
!    dr1dT = -Ufon * DtauDt
!
!    Mr = M * r1 / Rn
!!    dMrDx = M * (1d0 + Ufon*dTauDx - dRndX**2) / Rn
!!    dMrDy = M * (Ufon*dTauDy - r1*r2 / Rn**2) / Rn
!!    dMrDt = -c*M/Rn * (M*dTauDt + r1**2*(1-dTauDt)/Rn) ! c*(1-dTauDt) ! Ufon * M * (1d0 - r1 / Rn**2) / Rn
!    dMrDx = M * (1d0 - (r1/Rn)**2) * dr1dX / Rn
!    dMrDy = M * (dr1dY - r1/Rn**2*(r1*dr1dY + r2)) / Rn
!    dMrDt = M * (1d0 - (r1/Rn)**2) * dr1dT / Rn
!
!    dMrDx = M * (dr1dX - r1/Rn*dRndX) / Rn
!    dMrDy = M * (dr1dY - r1/Rn*dRndY) / Rn
!    dMrDt = M * (dr1dT - r1/Rn*dRndT) / Rn
!
!    s = k * Rn
!    omegaT = omega * Tau
!    
!    cosT = cos(omegaT)
!    sinT = sin(omegaT)
!
!    Y0 = DBesY0(s)
!    J0 = DBesJ0(s)
!
!    dJ0ds = -DBesJn(1,s)
!    dY0ds = -DBesYn(1,s)
!    
!!    HQ = Y0 * sinT + J0 * cosT
!!    dHdsQ = dY0ds * sinT + dJ0ds * cosT
!!    HdQdTau = omega * (Y0 * cosT - J0 * sinT)
!
!    HQ = Y0 * cosT - J0 * sinT
!    dHdrQ = k * (dY0ds * cosT - dJ0ds * sinT)
!    HdQdTau = -omega * (Y0 * sinT + J0 * cosT)
!    
!!    p = A * (dHdrQ * dRndTau + HdQdTau - HQ * dMrDtau / Mr1) * dTauDt / Mr1
!    p = A * (dHdrQ * dRnDt + HdQdTau * dTauDt + HQ * dMrDt * dTauDt) * dTauDt
!    
!    p = A * (HQ) * dTauDt
!    
!    U(1) = -A * (dHdrQ * dRndX + HdQdTau*dTauDx  + HQ*dMrDx/Mr1) / (Mr1*RoFon) + Ufon
!    U(2) = -A * (dHdrQ * dRndY + HdQdTau*dTauDy  + HQ*dMrDy/Mr1) / (Mr1*RoFon)
!
!    Rho = p / c**2 + RoFON
!
!    p = p + PFON

!    return
!    
!!    M = UFON / c
!!    Mr = UFON/c * r(1) / rMod
!!    beta = sqrt(1d0-M**2)
!!!    k = omega * sqrt(1d0/c**2 + 1d0/UFON**2)
!!    k = omega/c ! *beta) ! * sqrt(UFON**2-1d0) / UFON ! sqrt((omega)**2-1d0/UFON**2)/c
!!    dr = 1d-5
!!    r(1) = r(1) / beta ! - UFON
!!    time = beta * time + r(1)*M/c ! *beta)
!
!    r(1) = r(1) + Ufon * time
!    k = omega/c
!
!    rMod = sqrt(sum(r*r))
!    
!    omegaT = omega * time / beta !  - r(1)/UFON)
!    s = k * rMod !/ (1d0 + Mr)
!    
!    BesY0s = DBesY0(s)
!    BesJ0s = DBesJ0(s)
!    
!!    p = A*(cos(omegaT)*BesY0s + sin(omegaT)*BesJ0s)
!    p = A*(cos(omegaT)*BesY0s - sin(omegaT)*BesJ0s)
!    
!!    pt = A*(sin(omegaT)*BesY0s - cos(omegaT)*BesJ0s) * omega
!!    
!!    ptt = A*(cos(omegaT)*BesY0s + sin(omegaT)*BesJ0s) * omega**2
!    
!!    p = A*(cos(omegaT)*BesY0s - sin(omegaT)*BesJ0s)
!    
!!    p = 0.25d0*A*(sin(omegaT)*BesY0s + cos(omegaT)*BesJ0s)
!!    p = -0.25d0*A*(sin(omegaT)*BesY0s + cos(omegaT)*BesJ0s)
!
!!    do i = 1,2
!!        r1 = r
!!        r1(i) = r1(i) + dr
!!        rMod1 = sqrt(sum(r1*r1))
!!        drMod = rMod1 - rMod
!!        s2 = k * rMod1
!!        gradBesY0(i) = (DBesY0(s2)-BesY0s) / dr
!!        gradBesJ0(i) = (DBesJ0(s2)-BesJ0s) / dr
!!    end do
!
!    BesJ1s = DBesJn(1,s)
!    BesY1s = DBesYn(1,s)
!
!!    gradBesY0(1) = - k * r(1) * BesY1s / rMod
!!    gradBesY0(2) = - k * r(2) * BesY1s / rMod
!!    
!!    gradBesJ0(1) = - k * r(1) * BesJ1s / rMod
!!    gradBesJ0(2) = - k * r(2) * BesJ1s / rMod
!
!    gradBesY0(1) = - k * r(1) * BesY1s / rMod
!    gradBesY0(2) = - k * r(2) * BesY1s / rMod
!    
!    gradBesJ0(1) = - k * r(1) * BesJ1s / rMod
!    gradBesJ0(2) = - k * r(2) * BesJ1s / rMod
!    
!!    p = RoFon * A / beta * (omega/beta *(sin(omegaT)*BesY0s - cos(omegaT)*BesJ0s) - &
!!        UFON * (gradBesY0(1)*cos(omegaT) - gradBesJ0(1)*sin(omegaT)) )
!!        
!!    p = -RoFon * A / beta *(omega * BesJ0s*cos(omegaT) - UFON * k * r(1) * BesY1s / rMod * sin(omegaT))
!    
!    ! p == phi
!    Phi = A * ( BesY0s * cos(omegaT) - BesJ0s * sin(omegaT) )
!    
!    PhiS = - A * ( BesY1s * cos(omegaT) - BesJ1s * sin(omegaT) )
!    
!    PhiX = PhiS * k * r(1) / rMod
!    PhiY = PhiS * k * r(2) / rMod
!    
!    PhiT = - omega * A * ( BesY0s * sin(omegaT) + BesJ0s * cos(omegaT) )
!    
!!    p = - RoFon * (PhiT - UFON * PhiX) / beta
!    p = - RoFon * PhiT
!    
!!    U(1) = (M * PhiT / c + PhiX) / beta
!!    U(2) = PhiY
!
!    U(1) = PhiX + Ufon
!    U(2) = PhiY
!
!
!!    gradBesY0(1) = -BesY1s / rMod * (k * r(1) - (UFON / c - Mr*r(1) / rMod ) / (1d0 + Mr)**2)
!!    gradBesY0(2) = -BesY1s / rMod * (k * r(2) - (         - Mr*r(2) / rMod ) / (1d0 + Mr)**2) 
!!    
!!    gradBesJ0(1) = -BesJ1s / rMod * (k * r(1) - (UFON / c - Mr*r(1) / rMod ) / (1d0 + Mr)**2)
!!    gradBesJ0(2) = -BesJ1s / rMod * (k * r(2) - (         - Mr*r(2) / rMod ) / (1d0 + Mr)**2) 
!    
!    gradBesY0(3) = 0d0
!    gradBesJ0(3) = 0d0
!!    U(1:3) = - A*(cos(omegaT)*gradBesY0(1:3) + sin(omegaT)*gradBesJ0(1:3) ) / (ROFON*4*omega)
!
!!    rMod1 = rMod + dr
!!    s = k * rMod1
!!    gradBesY0 = (DBesY0(s)-BesY0s) / dr
!!    gradBesJ0 = (DBesJ0(s)-BesJ0s) / dr
!
!    do i = 1,3
!!        U(i) = A*(-sin(omegaT)*gradBesY0(i) + cos(omegaT)*gradBesJ0(i) ) / (ROFON * omega * beta)
!!        U(i) = A*(sin(omegaT)*gradBesY0(i) + cos(omegaT)*gradBesJ0(i) ) / (ROFON * c * beta) ! * beta)
!
!    !    U(i) = -r(i)*A*(-sin(omegaT)*gradBesY0(i) + cos(omegaT)*gradBesJ0(i) ) / (ROFON* rMod * omega)
!!        U(i) = r(i)*A*(-cos(omegaT)*BesY1s + sin(omegaT)*BesJ1s ) / (ROFON* c * rMod)
!
!!        Ut(i) = -A*(cos(omegaT)*gradBesY0(i) + sin(omegaT)*gradBesJ0(i) ) / (ROFON)
!!        Ut(i) = r(i)*A*(cos(omegaT)*gradBesY0(i) + sin(omegaT)*gradBesJ0(i) ) / (ROFON* rMod)
!!        Ut(i) = omega * r(i)*A*(sin(omegaT)*BesY1s + cos(omegaT)*BesJ1s ) / (ROFON* c * rMod)
!!        U(i) = r(i)*A*(cos(omegaT)*BesY1s + sin(omegaT)*BesJ1s ) / (ROFON* c * rMod)
!
!!        U(i) = r(i)*A*(cos(omegaT)*gradBesY0 - sin(omegaT)*gradBesJ0 ) / (ROFON*4*omega * rMod)
!!         U(i) = r(i)*A*(cos(omegaT)*BesY0s + sin(omegaT)*BesJ0s) / (ROFON*c*rMod)
!!         U(i) = r(i)*A*(cos(omegaT)*gradBesJ0 + sin(omegaT)*gradBesY0 ) / (ROFON*4*omega * rMod)
!!        U(i) = r(i)*A*(cos(omegaT)*gradBesY0 - sin(omegaT)*gradBesJ0 ) / (ROFON*4*omega * rMod)
!!        U(i) = r(i)*A*(cos(omegaT)*gradBesY0 - sin(omegaT)*gradBesJ0 - gradBesY0) / (ROFON*4*omega * rMod)
!!        U(i) = r(i)*A*(cos(omegaT)*gradBesY0 - sin(omegaT)*gradBesJ0 ) / (ROFON*4*omega * rMod)
!    end do
!    
!!    U(1) = U(1) * (UFON / (c*rMod) - UFON*r(1)**2 / (c*rMod**3) ) / (1d0 + Mr)**2
!!    U(2) = U(2) * (- UFON*r(1)*r(2) / (c*rMod**3) ) / (1d0 + Mr)**2
!    
!!    U(1) = UFON + (U(1) - omega*M/(c*beta) * (sin(omegaT)*BesY0s + cos(omegaT)*BesJ0s)) / beta ! / beta
!!
!!    U(1) = UFON + A/beta * (-k * r(1) * BesY1s / rMod * sin(omegaT) + omega * M/c * BesJ0s*cos(omegaT))
!!    U(2) = -A * k * r(2) * BesY1s / rMod * sin(omegaT)
!
!!    U(1) = (M/c * PhiT + PhiX) / beta + UFON
!!    U(2) = PhiY
!
!    U(1) = PhiX + Ufon
!    U(2) = PhiY
!
!
!    Rho = p / c**2 + RoFON
!
!    p = p + PFON
!    
!    
!!!    rMod1 = rMod - dr
!!!    s = k * rMod1
!!!    gradBesY01 = (BesY0s-DBesY0(s)) / dr
!!!    gradBesJ01 = (BesJ0s-DBesJ0(s)) / dr
!!!
!!!    ddBesY0 = (gradBesY0-gradBesY01) / dr
!!!    ddBesJ0 = (gradBesJ0-gradBesJ01) / dr
!!    
!!    do i = 1,2
!!        r1 = r
!!        r1(i) = r1(i) + dr
!!        rMod1 = sqrt(sum(r1*r1))
!!        s2 = k * rMod1
!!
!!        r1 = r
!!        r1(i) = r1(i) - dr
!!        rMod1 = sqrt(sum(r1*r1))
!!        s1 = k * rMod1
!!
!!        ddBesY0(i) = ((DBesY0(s2)-BesY0s)/dr-(BesY0s-DBesY0(s1))/dr) / dr ! (DBesY0(s2)-2d0*BesY0s+DBesY0(s1)) / dr**2
!!        ddBesJ0(i) = ((DBesJ0(s2)-BesJ0s)/dr-(BesJ0s-DBesJ0(s1))/dr) / dr ! (DBesJ0(s2)-2d0*BesJ0s+DBesJ0(s1)) / dr**2
!!!        gradBesY0(i) = (DBesY0(s)-BesY0s) / dr
!!!        gradBesJ0(i) = (DBesJ0(s)-BesJ0s) / dr
!!    end do
!!    ddBesY0(3) = 0d0
!!    ddBesJ0(3) = 0d0
!!    
!!    gradP(1:3) = -A*(-cos(omegaT)*gradBesY0(1:3) - sin(omegaT)*gradBesJ0(1:3))
!!    divU = -A*(-sin(omegaT)*sum(ddBesY0(1:3)) + cos(omegaT)*sum(ddBesJ0(1:3)) ) / (ROFON* omega)
!!    
!!    !p       = A*(-cos(omegaT)*BesY0s -            sin(omegaT)*BesJ0s)
!!    laplaceP = A*(-cos(omegaT)*sum(ddBesY0(1:3)) - sin(omegaT)*sum(ddBesJ0(1:3)))
!!    
!!!    divU = -2*A*(-sin(omegaT)*gradBesY0 + cos(omegaT)*gradBesJ0 ) / (ROFON* rMod * omega) - &
!!!        A*(-sin(omegaT)*ddBesY0 + cos(omegaT)*ddBesJ0 ) / (ROFON* omega)
!!    
!!    Err(1:3) = Ut(1:3) * RoFon - gradP(1:3)
!!    
!!    ErrPt = Pt / c**2 + RoFon * divU
!!    
!!    ErrPt1 = Pt + 1.4d0 * PFon * divU
!!    
!!    ErrPtt = Ptt/c**2 - laplaceP
!!
!!    dJ1 = -BesJ0s + DBesJn(1,s) / s ! -0.5*(BesJ0s+DBesJn(2,s))
!!    dY1 = -BesY0s + DBesYn(1,s) / s ! -0.5*(BesY0s+DBesYn(2,s))
!!
!!    s1 = s - dr
!!    s2 = s + dr
!!
!!    dJ1c = ((DBesJ0(s2)-BesJ0s) / dr - (BesJ0s - DBesJ0(s1))/ dr ) / dr
!!    dY1c = ((DBesY0(s2)-BesY0s) / dr - (BesY0s - DBesY0(s1))/ dr ) / dr
!!!    dY1c = (DBesY0(s1)-2d0*BesY0s+DBesY0(s2)) / dr**2
!!    
!!    dJ0s = (DBesJ0(s2)-BesJ0s) / dr
!!    dY0s = (DBesY0(s2)-BesY0s) / dr
!!    
!!    dJ0s1 = (BesJ0s - DBesJ0(s1))/ dr
!!    dY0s1 = (BesY0s - DBesY0(s1))/ dr

end subroutine CalcQ2D_Hankel


! Analytic 3D solution
subroutine CalcQ3D_ConvectedMonopole(r_in,time_in,p,U,rho)
use Data_cfoam, only: PFon,RoFon, UFON
implicit none
    real(8),intent(in) :: r_in(3),time_in
    real(8),intent(out) :: p,U(3), rho

    integer :: i
    real(8) :: A, omega, c, k, s, omegaT, dr, rMod, rMod1, drMod !, r1(3)
    real(8) :: BesY0s, BesJ0s, r(3)
!    real(8) :: gradBesY0(3),gradBesJ0(3)
    real(8) :: gradBesY0(3),gradBesJ0(3), BesJ1s, BesY1s, ddBesY0(3),ddBesJ0(3)
    real(8) :: Ut(3), gradP(3), Err(3), Pt, divU
    real(8) :: gradBesY01, gradBesJ01 !, ddBesY0, ddBesJ0
    real(8) :: s1, ErrPt,ErrPtt,Ptt,laplaceP, ErrPt1
    real(8) :: Pxx,Pyy,Pzz, dY1,dJ1,dY1c,dJ1c, s2, dJ0s,dY0s, dJ0s1,dY0s1
    real(8) :: time, M, beta
    real(8) :: Phi, PhiS,PhiT,PhiX,PhiY,PhiZ
    
    real(8) :: beta2, Mx1c
    
    real(8) :: Tau,dTauDt,dTauDx, dTauDy, D
    real(8) :: r1,r2, Rn, dRnDt, dRndX, dRndY
    real(8) :: J0,Y0, dJ0Dr,dY0Dr, cosT,sinT
    real(8) :: Mr, dMrDt, dMrDx, dMrDy, Mr1
    
    real(8) :: HQ, dHdrQ, HdQdTau
    real(8) :: dJ0ds,dY0ds
    
    real(8) :: dr1dX, dr1dy, dr1dT    
    
    real(8) :: r1_,r2_,Rn_,Mr_,dMrDtau_
    real(8) :: TauT,RnT,r1T, TauX,RnX,r1X
    real(8) :: Pi

    Pi = 3.1415926535897932384626433832795

    A = 200d0
    omega = 2d0*Pi ! 1d0

    r = r_in
    time = time_in

    c = sqrt(1.4*PFon/RoFon)
    M = Ufon/c

    beta2 = 1d0-M**2    
    beta = sqrt(beta2)
!    Rn = sqrt((r(1))**2 + (r(2)**2+ r(3)**2)*beta2)
    Rn = sqrt((r(1)/beta)**2 + r(2)**2+ r(3)**2)
!    tau = omega * (time + (M*r(1) - Rn)/(c*beta2))
    tau = omega * (time + (M*r(1)/beta - Rn)/(c*beta))

    Phi = A * sin(tau) / (beta*Rn)
    PhiT = A * omega* cos(tau) / (beta*Rn)
    
    PhiX = A * ( omega*cos(tau) * (M - r(1)/(beta*Rn))/c - sin(tau) * r(1) / (beta*Rn**2) ) / (beta**2*Rn)
    PhiY = -A * ( omega*cos(tau)/(beta*c) + sin(tau)/Rn ) * r(2) / (beta*Rn**2)
    PhiZ = -A * ( omega*cos(tau)/(beta*c) + sin(tau)/Rn ) * r(3) / (beta*Rn**2)
    
    
!    Phi = A * cos(tau) / Rn
!    PhiT = -A * omega* sin(tau) / Rn
!
!    PhiX = A * ( omega*sin(tau) * (-M + r(1)/Rn)/(c*beta2) - cos(tau) * r(1) / Rn**2 ) / Rn
!    PhiY = A * ( omega*sin(tau)/(c*beta2) - cos(tau)/Rn ) * r(2) / (beta2*Rn**2)
!    PhiZ = A * ( omega*sin(tau)/(c*beta2) - cos(tau)/Rn ) * r(3) / (beta2*Rn**2)
!    omega*beta * cos(tau) * (M/(c*beta) - r(1)/(c*beta2*Rn)) - r(1) * sin(tau) / (beta*beta2*Rn**3)

    p = -RoFon * (PhiT + Ufon * PhiX)
    
    U(1) = PhiX + Ufon
    U(2) = PhiY
    U(3) = PhiZ

    Rho = p / c**2 + RoFON

    p = p + PFON

    return
    
!contains
!
!    subroutine CalcTauRn(r,time,Tau,Rn,r1)
!    implicit none
!        real(8),intent(in) :: r(3),time
!        real(8),intent(out) :: Tau,Rn,r1
!        
!        real(8) :: D
!        
!        D = (r(1) - Ufon*time)**2 + beta2*r(2)**2
!        D = sqrt(D)
!        
!        Tau = (c*time - M*r(1) - D) / (c*beta2)
!        Rn = c*(time-Tau)
!        r1 = r(1)-Ufon*Tau
!    end subroutine CalcTauRn

end subroutine CalcQ3D_ConvectedMonopole

subroutine BlowoffInitialCondition(r,p,U,Rho)
implicit none
    real(8),intent(in) :: r(3)
    real(8),intent(out) :: p,U(3),Rho

    real(8) :: x, x1,x2,x3
    real(8) :: p1,u1,Rho1, p2,u2,Rho2, p3,u3,Rho3, p4,u4,Rho4

    x1 = 0.35d2
    x2 = 0.5d2
    x3 = 0.55d2

    p1=1d2
    u1=0d0
    Rho1=5d-1

    p2=2d2
    u2=-5d-1
    Rho2=7.5d-1

    p3=4d2
    u3=5d-1
    Rho3=1d0

    p4=1d2
    u4=0d0
    Rho4=5d-1

    x = r(1)

    if(x.lt.x1)then
        p = p1
        U(1) = u1
        Rho = Rho1
        return
    end if

    if(x.lt.x2)then
        p = p2
        U(1) = u2
        Rho = Rho2
        return
    end if

    if(x.lt.x3)then
        p = p3
        U(1) = u3
        Rho = Rho3
        return
    end if

    p = p4
    U(1) = u4
    Rho = Rho4
    return

end subroutine BlowoffInitialCondition
